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Abstract 



A quantum fluid dynamic(QFD) control formulation is presented for optimally 
manipulating atomic and molecular systems. In QFD the control quantum system is ex- 
pressed in terms of the probability density p and the quantum current j. This choice of 
variables is motivated by the generally expected slowly varying spatial-temporal depen- 
dence of the fluid dynamical variables. The QFD approach is illustrated for manipulation 
of the ground electronic state dynamics of HCl induced by an external electric field. 

I. INTRODUCTION 

Manipulating the outcome of quantum dynamics phenomena by a properly tailored 
external control is a topic of increasing activity [1-15]. Problems where the external con- 
trol is electromagnetic have received the most attention, although other applications arise 
as well. Various implementations of quantum control have been experimentally realized 
[16-25]. A variety of control strategies have been suggested, and optimal control thc- 
ory(OCT) provides the most general framework for acheiving field designs. Such designs 
will generally require further refinement in the laboratory through learning techniques to 
overcome design uncertainities [17-20] . A basic difficulty in attaining the control designs 
is the computational effort called for in solving the time-dependent Schroedinger equa- 
tion, often repeatedly in an iterative fashion. This paper introduces the quantum fluid 
dynamic(QFD) control formulation to simplify this task. Calculations have shown [28, 
29] that QFD is capable of being much more efficient than conventional methods(e.g., 
FFT propagation), and this savings should carry over to the control design task. This 
paper will show how OCT can be combined with QFD. 

The theoretical basis for dynamic control [2-13] is to create non-stationary states of 
one's choice, by optimally designing the control field. Typically, the problem is posed as 
seeking an optimal field to drive a quantum wave packet to a desired target at a chosen 



2 



time t=T. In the traditional approach [3-10] to quantum optimal control an objective 
design functional J is defined, which depends on the system wave function, a wave function 
like Lagrange multiplier, and the external field. Minimization of the objective functional 
leads to the identification of external field(s) capable of delivering a specific outcome. 
This process requires solving for the complex oscillatory wave function and the similarly 
behaved Lagrange multiplier, and due care is needed for their proper representation in a 
suitable basis for capturing their behaviour. Often the rapidly varying spatio-temporal 
behaviour of these functions necesitates the use of many unknowns in the basis. 

This paper explores an alternative formulation for OCT to design the electric field. 
The formulation is based on the fluid dynamic view point of quantum mechanics [26- 
29], which by-passes the typically oscillatory nature of the wave function to exploit the 
generally smooth behaviour of the real density and the quantum current variables. Recent 
illustrations have demonstrated the smooth spatial and temporal nature of the variables 
and the ability to discretize them on a relatively small number of grid points [28, 29]. 
As background in section 2 we give a brief summary of the QFD formulation of quantum 
mechanics. Section 3 presents OCT within the framework of QFD for designing an electric 
held to meet a speciflc objective. Section 4 apphes the OCT-QFD formulation for the 
manipulation of HCl. Section 5 concludes the paper. 



II. QUANTUM FLUID DYNAMICS 

The treatment below considers a single particle of reduced mass m, but 
the QFD formulation has an immediate extension to many particles. The time-dependent 
Schrodinger equation is given by 

[-|^V^ + V + K.*]^(x, t) = ^^^*(x, t) (1) 

where V typically confines the particle in a locale and Vg^t is the control taken here as 
— /i(x) • E{t) with E(t) being the electric field and /i(x) the dipole moment. Substituting 
^(x, t) = A(x,t)e^^^^'^^^^, where A and S are real functions, into Eq.(l) and separating 
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the real and imaginary parts one easily obtains two equations. The imaginary part yields 
the continuity equation 

|^ + V-(pv) = (2) 
and the real part the following equation for the phase S 

where = ^ + + K with V, = = -^[V^/npV^ + {Vlnp^^n P = 

and v=— . 

m 

Equation (3) has the form of the classical Hamilton-Jacobi equation with an extra 
'quantum potential' term Vg. This equation can be transformed into one for the evolution 
of the velocity vector v by taking the gradient to give 

^v = -(vV)v-lv(K//) (4) 

Defining the quantum current as 

j(x,i) = --/m[**(x,i)V*(x,i)] =p(x,t)v(x,i), 

one readily obtains the equation of motion for j by substitution of || and from Eqs.(2) 
and (4) as 

|j = -v(V-j)-(j-V)v-^VK// (5) 

Eqs.(2) and (3) or (2) and (5) describe the motion of a quantum particle within the QFD 
formulation of quantum mechanics. The motion of a quantum particle is governed by 
the current vector j and the density p in Eqs.(2) and (5). Although the QFD equations 
resemble those of classical fluid dynamics, their quantum identity prevails due to the 
presence of the potential Vg which has no classical analogue. Equivalently, the QFD 
equations may be viewed as those of a "classical" fluid with a highly non- linear constitutive 
law prescribed by K//. Various Eulerian or Lagrangian means can be exploited to solve 
the QFD equations [28,29], and available fluid dynamics codes may be adopted to treat 
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these equations [30]. The essential simphfying feature of the QFD equations is that p 
and j or p and S are often slowly varying, which is evident from quantum dynamics 

calculations [28,29], thereby permitting relatively coarse gridding. Despite the non-linear 
nature of the QFD equations, the general smoothness of p and j observed lead to significant 
computational savings [28,29]. 

III. CONTROL EXPRESSED WITHIN QUANTUM FLUID 

DYNAMICS 

Quantum OCT seeks the design of an external field to fulfill a particular dynamical 
objective. This section will provide the working equations for OCT-QFD to design an 
optimal electric field that drives a quantum wave packet to a desired objective at the 
target time t=T. The OCT-QFD formulation could be expressed in the usual way in 
terms of the Schrodinger equation where QFD would only act as a solution procedure. 
Here we will present a general approach by writing OCT directly in terms of QFD. As 
an example the control of a non-rotating diatomic molecule will be used as a simple 
illustration of the concepts. The treatment of a fully general target expectation value 
Qj, =< *(r)|©|^'(r) > may be considered with QFD, but here we will only treat the 
common case where the operator ©(x) is only position x dependent. Then the goal is to 
steer Oy 

©r= r Q{x)p{x,T)dx (6) 

Jxi 

as close as possible to the desired value Q"^. 

The active spatial control interval is taken os xi < x < Xr over the time < i < T that 
the control process occurs. We desire to minimize the cost functional Jcost = Jtarget + J field 
where Jtarget and Jfieu are given by 

Jtarget = -;<^x{Qt - Q'^f and J field = -UJe E'^{t)dt (7) 

with cUe and cu^ being the positive weights balancing the significance of the two terms. 
The second term represents the penalty due to the fiuence of the external field. The 
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minimization of Jcost with respect to E(t) must be subject to the satisfaction of the 
equations of motion for p and j in Eqs.(2) and (5). We may fulfill this constraint by 
introducing the unconstrained cost functional as 

7 T l'^ . , ^^dp(x,t) dj(x,t), , , 

where \i{x,t) and \2{x,t) are Lagrange's multiplier functions. 

An optimal solution satisfies 5J = 0, which is assured by setting each of the functional 
derivatives with respect to Ai, A2, p, j and E to zero. The first two, i.e., the functional 
derivatives with respect to Ai and A2 regenerate the QFD equations in Eq.(2) and (5). 
The three others are obtained in the forms : 

^ + ^{X2Vx) + S,[p,j,X2]^0 (9) 

' + ^(Alt^A) - \2^{V + V,{\2) + Ve,t) + ^2[P, J, A2] = (10) 



and 



dt dx dx 



/ X2{x,t)p{x,t)—p{x)dx + u;eE{t)^0 (11) 

Jxi OX 



5E{t) Jxi ' ' ' 'dx 



where 



S.=2l^ (12) 

p OX 



[tT72^(M] 



Aim?p^/'^ dx^ p^/^ dx 
1/2 d_, 

+ ^^2p3/2Q^2p Qx^^^P' 



and 
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2m A^/' 



2m 



(14) 



The corresponding final conditions are 



c^.[©T-e'^]e(x)-Ai(a;,r) = 



(15) 



and 



A2(x,r) = 



(16) 



Several other constraint expressions can be obtained by using equivalent forms of the 
continuity and dynamical equations. The form presented above is used in the subsequent 
numerical calculations. An alternative form in multi-dimensions symmetric between the 
QFD and Lagrange multipher functions is presented in the Appendix. 

The equations (9) and (10) for A2 and Ai respectively ressemble that of p and j with 
the only difference being the extra source terms Si and 82- The source terms depend on p 
and j. v\ in the above equations is the 'velocity' associated with the Lagrange's multiplier 
and is given as v\ = There are now two different quantum potential terms, one of 
which is a function of p{x,t) and the other is a function of X2{x,t). In this formaUsm 
the evolution of Xi{x,t) takes place by Vg{X2) as well as the difference of the two types 
of quantum potential. In obtaining the above equations we have standardly assumed no 
variation of either p{x,0) or j(x,0). Thus, we start from the initial value of p{x,0) and 
j(x,0) to solve Eqs.(2) and (5). Eqs.(9) and (10) can be solved for X2{x,t) and Xi{x,t) 
by integrating backward from time T using Xi{x,T) and X2{x,T) given in Eqs.(15) and 
(16) respectively. The equations (2), (5), (9) and (10) are non-linear thereby calling for 
iteration to solve(cf., the algorithm in Section 4). Finally the desired control electric field 
is given from Eq.(ll) as 




(17) 
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4 APPLICATION TO HCL 



The OCT-QFD formulation will be applied to manipulating the vibrational motion of 
HCl on the ground electronic state. The initial density p{x,0) = {^{x)]"^ was obtained 
from solving for the vibrational state from the equation 

^ + V{x)^{x) = E^{x) (18) 

m ox^ 

using the Fourier grid Hamiltonian method [31,32] where m is the reduced mass of the 
HCl molecule and V{x) is the truncated polynomial presented by Olgilvie [33] 



[V{x) 



«i(2f^f;)'[l + EL«.(2f=f;^V&i for x<4 

A[l - tanh{x - A)f''^ for A<x< 6.5 ] (19) 
for X > 6.5 



where Xe=2.4086 a.u. is the equilibrium bond length of HCl. The parameters in a.u. 
entering the potential function are ai = 0.961914, 02 = —1.362999, 03 = 0.86675, 04 = 
-0.49804, as = 0.1727, ae, = 0.2687, 0,7 = -1.977, ag = 2.78, ag = 4.89, h = 0.169695 
and A — —4.85 x 10~^. Since 'i/{x) is a stationary real function we have zero initial flux 
j(x,0)=0. The initial p{x, 0) is nearly a Gaussian packet centered around Xg. The dipole 
function for HCl is given by [34] 

li{x) = ci[g{x) + Cig^x) + csg^ix)] (20) 

where g{x) — 1 — tanh{f3{x — xa)) and the parameters are ci = 0.279, C2 = —0.905, 
C3 = 1.029, /3 = 0.687 and xa — 2.555. The following steps were carried out for imple- 
mentation of the present OCT-QFD algorithm : 



Step 1: Make an initial guess for the electric Geld E(t), which was zero in the 

present calculations. 

Step 2: Solve the coupled equations, viz., Eq.(2) and (5) for p{x,t) and 
j(x,t) respectively starting from p{x, 0) and j(x,0). The solution was 
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achieved here by using the Flux- corrected transport (FCT) algorithm 
[35] modihed for the purpose of solving the QFD equations [28]. In 
doing so, we adopt the Eulerian numerical scheme. 

Step 3: Evaluate the final value for Ai(x, T) given by Eq.(15) and set 
X^{x,T)=0 byEq.(16). 

Step 4: Solve Eqs.(9) and (10) for \2{x,t) and Xi{x,t), respectively, by back- 
ward propagation using the same method as in step 2. Equations 
(9) and (10) have source terms which depend on p{x,t) and j(x,t) 
calculated from step 2. 

Step 5: Calculate the difference between the left and right sides of Eq.(16) 
for use in the conjugate gradient method [36] and calculate Jcost from 
Eq.(7). 

Step 6: Iterate steps 2 to step 6 until acceptable convergence is met. 

The spatial range of the calculation was < a; < 12 a.u., and the time interval was 
< t < T with T=2000 a.u. The total number of spatial mesh points is 64 which gives 
Ax = 0.1875 a.u. Similarly, the total number of time steps was 2048, which corresponds 
to At = 0.9765 a.u. No special effort was made to optimize the grid points, as the purpose 
here is to demonstrate the QFD-OCT formulation. The weight Ue in Eq.(7) was taken as 
|, and = 1000. The target operator was Q = x and 6*^ = S.Oa.u.. 

Figure |T] shows the control field in atomic units. The slightly non-zero values of the 
field at the beginning and end could be arrested by placing additional costs if desired. 
This pulse excites several vibrational states(not shown here) mainly by a sequence of 
single quantum transitions. Figure |^ shows the average distance < x > as a function of 
time. The desired control value of < x >=3.0 a.u. at T is obtained through oscillatory 
motion of the packet. The packet is distorted in shape (not shown) while approximately 
retaining its original variance during the evolution. During the optimization process the 
total integrated probability density remained at unity up to a deviation of 10~^. The 
iteration algorithm takes 10 steps to achieve the results shown here at 2 CPU mins. on 
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an IRIX Silicon Graphics Machine (Release 6.1). Within numerical precision the results 
were the same as obtained by solving the original Schrodinger equation. 

5. CONCLUSION 

This paper presents a new QFD based approach for carrying out the optimal design of 
control fields with an illustration for the maniputation of the HCl molecule. Our previous 
work [28] shows the typical smooth and monotonic behaviour of the fluid dynamical 
variables, viz., S and v as opposed to the typical oscillations in the wave functions where 
the hamiltonian was time independent. In the present case where the system is driven 
with an optimal time-dependent external field we have calculated the spatial dependence 
of j, p, S and \1/ at t=T shown in Fig.^. The fluid dynamical variables(Fig.0 curves (a), 
(b) and (c)) used in the present method are relatively slowly varying spatial functions 
compared to the wave function(Fig.|^, curve(d)) which apparently enhances the efficiency 
and the numerical saving of the present approach to controlling dynamics. 

Although the illustration was for one dimension the QFD technique is directly extend- 
able to higher dimensions, and a QFD wave packet calculation in four dimension has 
already been performed [28]. The alternating direction method can effectively be used 
with QFD for high dimensions. Comparison with FFT propagation has been performed 
for two dimensional systems [29] , showing that QFD is capable of providing a considerable 
increase in efficiency(i.e., by a factor of 10 or more). Regardless of the dimension, the key 
advantage of OCT-QFD arises from the expected smooth nature of QFD variables. A spe- 
cial circumstance will arise if the control "exactly" leads to a bound state with nodes that 
fully separates one spatial region from another. In practice placing a lower limit on the 
density of the order of the machine precision overcomes such difficulties. Future studies 
need to explore the full capabilities of the computational savings afforded by OCT-QFD. 
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APPENDIX 

Two additional forms for the cost functional and associated initial/final 

conditions 

The forms here are presented for reference as an alternative QFD approach. They have 
the advantages of simplicity and of giving equations for the Lagrange multiplier in the 
same form as the dynamical equations. The formalism for deriving the quations is through 
the Euler equations corresponding to the minimization of 

1= [ r F{fJ„Vf,V'f)dtdV (Al) 
Jv Jt=o 

Here V denotes the volume in coordinate space. The corresponding Euler equations and 
conditions on time and space are 

df dt ^ ^dVf^^ ^dv'f-^ ^^-^^ 

Initial condition: /(x, 0) = /o(x); Final condition: = /t(x) (A. 3) 

Boundary conditions on dV: 

/(x, t) = /b(x, or n . - V • (^)] = 0; (A.4) 



n-V(/)=^B(x,i)or^ = (A.5) 



Lav/ ^ \dvf^ 

F 

Starting with the continuity and energy conservation equations given in Eqs.(2) and (3) 
in the text, we rewrite them as 

^, + V^.Z^ + ^ = (A6) 
m 2m 

ASt + A^^^ + VA- ^V'A = (A.7) 
2m 2m ^ ^ 

The use of the dynamical equations above in the cost functional in Eq.(8) becomes 
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J = \oJ^A^T - Q'f + r E\t)dt 
I Jt=0 



T V4 • V 9 
[[Ai(A + — + AV^S)/2m] {A.8) 

IV Jt=o m 

+[X2{ASt + A ^ +VA- iiE{t)A - i^V^A)]]dtdV 

where 0^ = /y 0(x)A^(x, T)(iK. The corresponding Euler equations are obtained from 
the formulas in {A2) for arbitrary variations of A, S, Ai, A2 and E{t) as 

At + v- VA = -AV-v/2 

5t + V • WS/2 = -V + iiEit) + — V'^M 

2m 

+ V • VAi = -[Ai^ + ^\^{V^A/A - V^A2)/A2)] (A.9) 
z Im 

A2t + V • VA2 = [A2^ + ^\^{y^A/A - V^AO/Ai)] 
z Im 



UeE{t) + I \2HAdV = 



Following the formulas given in Eq.(A3) to (A5), the corresponding initial and final con- 
ditions become 

>l(x,0) = Ao(x);5(x,0) = 5'o(x); 

Ai(x,r) + 2u;,(eT - e<^)^(x,r)e(x); A2(x,r) = o (a.io) 

The first two formulas in A. 9 are equivalent to the Schrodinger equation. They can 
be transformed into various QFD forms in terms of p, v and j as in Eqs.(2) to (5) in 
the main text. The third and fourth equations in A. 9 are the basic equations for the 
Lagrange multiplier functions. They are in the same flux conservation form as the QFD 
equations. Indeed, the third equation multiplied by Ai can be rearranged in the form of 
mass conservation for Ai = A^ as 
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Alt + V • (Aiv) = -[—XiX2{V^A/A - V'As/As)] 

The above derivation also can be obtained starting with the usual Schrodinger equation 
and its complex conjugate. Following this approach the cost functional below assures that 
the external field is real 

J = ^^.(Ot - Q''? + UJe r E\t)dt 

2 Jt=o 

- [ r [>^*['^^t + i-v* -v^- i^E{t)^] 

Jv Jt=Q 2m 
+A[-i** + ^V^** - V^* - iiE{t)^*]]dtdV 

With the substitution ^ = Aexp{iS), the cost functional reduces to the one in Eq.(A.8) 
with A = Ai + iX2- 
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t{a.u.) 



FIG. 2. The expectation value < x > shown as a function of time in a.u. The target value is 
< X >=3.0 a.u. at T=2000 a.u. 
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FIG. 3. Fluid dynamical variables, viz., p(x)(curve (a); dotted lines for the initial density 
and solid lines for the final density), j(x) (curve (b)), S(x) (curve (c)) shown as a function of x 
corresponding to t=T. Curve (d) shows the wave function(*(a;))(solid lines for real part and 
dotted lines for imaginary part) as a function of x corresponding to t=T. 
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